% Landau–Lifshitz方程
% 参考：W站“Landau–Lifshitz方程”词条
% 疑似有点太简单了
% Gitee Repo

clc
clear

gamma = 1; % 旋荷比，假定正电球
lambda = 0.01; % 阻力系数
dt = 0.01; % 时间步长

M0 = [1;0;1]; % 磁矩
M1 = M0;

B = [0;0;1]; %有效磁（感应）场强度

mode = 'sphere'; # 'sphere' or 'coil'

figure()

function h = SPHERE(x,y,z,r)
  [X,Y,Z] = sphere(20);
  h = surf(r*X+x,r*Y+y,r*Z+z);
  set(h,'EdgeColor','none')
  set(h,'CData',X)
end

function myrotate(obj, normvec)
    % 这段代码是ai写的，我不知道对不对...
    % 将图形对象 obj 的法向量从 (0,0,1) 旋转至目标法向量 normvec = [a, b, c]
    % 输入:
    %   obj: 图形对象句柄（如 surf, patch 等）
    %   normvec: 目标法向量 [a, b, c]（无需单位化）

    % 默认法向量
    current_normal = [0;0;1];

    % 计算旋转轴（叉积）
    rotation_axis = cross(current_normal, normvec);

    % 计算旋转角度（点积）
    rotation_angle = acos(dot(current_normal, normvec) / norm(normvec));

    % 处理特殊情况（normvec 平行或反平行于 (0,0,1)）
    if norm(rotation_axis) < eps
        if normvec(3) < 0 % 反平行 (0,0,-1)
            rotation_axis = [1;0;0]; % 任意轴（如 x 轴）
            rotation_angle = pi;       % 180度
        else % 平行 (0,0,1)，无需旋转
            return;
        end
    else
        % 归一化旋转轴
        rotation_axis = rotation_axis / norm(rotation_axis);
    end

    % 执行旋转（角度转换为度）
    rotate(obj, rotation_axis, rad2deg(rotation_angle));
    return
end

for tick = 1:5000
    M1 = M0 + dt*(gamma*cross(M0,B)-lambda*cross(M0,cross(M0,B)));
    M0 = M1;

    if mod(tick,50) == 0
        clf
        hold on
        axis equal
        axis([-2,2,-2,2,-2,2])

        quiver3(0,0,0,M1(1),M1(2),M1(3),0); %磁矩向量
        quiver3(0,0,0,B(1),B(2),B(3),0); % 磁（感应）场强度向量

        t = linspace(0,2*pi,30);

        if strcmp(mode,'coil')
            h = plot3(cos(t),sin(t),zeros(size(t)));
        elseif strcmp(mode,'sphere')
            h = SPHERE(0,0,0,1);
            rotate(h,[0;0;1],rad2deg(norm(0.5*M1*dt*tick)))
        else
            error('wrong mode')
        end
        myrotate(h,M0)

        xlabel('x')
        ylabel('y')
        zlabel('z')
        view(30,45)

        drawnow
        pause(0.01)
    end
end


